Transition from viscous to inertial regime in dense suspensions 
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Non-Brownian suspensions present a transition from Newtonian behaviour in the zero-shear limit 
to a shear thickening behaviour at a large shear rate, none of which is clearly understood so far. 
Here, we carry out numerical simulations of such an athermal dense suspension under shear, at 
an imposed confining pressure. This setup is conceptually identical to the recent experiments of 
Boyer and co-workers pQ. Varying the interstitial fluid viscosities, we recover the Newtonian and 
Bagnoldian regimes and show that they correspond to a dissipation dominated by viscous and 
contact forces respectively. We show that the two rheological regimes can be unified as a function 
of a single dimensionless number, by adding the contributions to the dissipation at a given volume 
fraction. 

PACS numbers: 83.80.Hj,47.57.Gc,47.57.Qk,82.70.Kj 



The rheology of amorphous materials such as emul- 
sions, foams, metallic glasses, suspensions or granular 
materials share a similar phenomenology close to the jam- 
ming transition at which viscosity diverges [IH1]- How- 
ever, the dynamics of these systems is not yet clearly 
understood and the establishment of a unified theory re- 
mains a challenging goal of out of equilibrium statistical 
physics. Following the pioneering work of Einstein [5], 
the common view on suspensions of particles in a fluid 
has long been to start from the dilute limit and to per- 
form an expansion in volume fraction </> [5J [TJ, with a 
particular emphasis on the effective interaction between 
particles mediated by the fluid. By contrast, recent stud- 
ies have started to view the rheology of dense suspensions 
from the other limit instead, in the framework of dense 
granular systems [2 IMT2]. The rheology of dense sus- 
pensions of solid particles in an isodense fluid of viscos- 
ity r)f is Newtonian at small shear rate 7 with a viscosity 
r/7 diverging as r]f(4> c — <j))~^ , as the particle volume 
fraction goes to its critical value <p c . The measured ex- 
ponent (3 ranges between 2 and 3 [TJ . Mean field 
theory assuming a dissipation dominated by lubrication 
films separating particles predicts an exponent j3 = 1 [S] . 
By contrast, numerical simulations assuming that dissi- 
pation is due to the nonaffine displacement of particles 
give the exponent /? ~ 2.2 [TUHH). There they relate the 
zero-shear viscosity, a macroscopic dynamical observable, 
to a microscopic observable: the variance of the nonaffine 
velocity /displacement [TO]. The latter is itself related to 
the geometry of the contact network [11] . 

While most fluids shear thin, it was first shown by 
Bagnold [16] that suspensions exhibit shear thickening 
when the volume fraction <f> is kept constant: their ap- 
parent viscosity increases with the shear rate. However, 
the conditions for such a property to emerge still remain 
controversial |12j . In particular, as recently emphasized 
[I] , suspensions exhibit shear thinning when the confining 
pressure P p is controlled and kept constant, a property 
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FIG. 1: Fraction of the power dissipated by contact forces 
(r cont ), viscous drag plus Archimedes forces (Fd+A), lubri- 
cation forces (Ti u br) and fluid viscosity (Ffl u id = f?/7 2 ), as a 
function of the ratio I / J for a fixed value of the volume frac- 
tion (<j> ~ 0.78). Solid lines are the best fits to the expression 
j+ai' 2 ' w ith Ci as fitting parameters, for the three last dissi- 
pation components. Insets: Friction coefficient /1 against I/J 
(solid curve illustrates the average p) and a schematic of the 
numerical setup. 



reminiscent of dry granular materials. 

In this Letter, we use discrete element simulations of 
non-Brownian particles interacting with a continuum vis- 
cous fluid to show that the rheology of suspensions at 
finite shear rate can be unified with the Newtonian qua- 
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FIG. 2: Friction coefficient of the suspension p — r/P p as a 
function of the particle volume fraction <j>, for different values 
of// J. Inset: Same as the main plot but for frictionless grains 
(/Hp = 0). The solid lines are the best fit by Eq. (??) for 
fi p — 0.4 and the dashed show the critical values for p p — 0.4. 



sistatic limit. More precisely, Boyer et al. pQ have re- 
cently shown that the rheology of a suspension approach- 
ing the zero-shear limit can be rewritten as a frictional 
law of the form r = pj(J)P p and <f> — <fij(J), where 
J = rifj/PP is the viscous number comparing viscous 
stresses to the confining pressure. In the inertial Bag- 
noldian regime, the flow is characterized by the inertial 
number / = W pj 2 d 2 /P p , with a subsequent rheology of 
the form r = pi(I)P p and <j> = 4>i(I). We show here that 
the contributions to the dissipation can be added at fixed 
(j), which results in a unique rheology r = fJ-(K) and 4>(K) 
controlled by the dimensionless number K = J + al 2 , 
where a is a constant of order 1 encoding the details of 
dissipative mechanisms. 

Numerical model. - We consider a two-dimensional sys- 
tem constituted of ~ 10 3 circular particles of mass 
and diameter a\, with a ±50 % polydispersity. The shear 
cell is composed by two rough walls, created by gluing 
together two dense layers of grains, with periodic bound- 
ary conditions along the direction x parallel to the walls. 
The position of the walls is controlled to ensure a con- 
stant normal stress P p and a constant mean shear ve- 
locity 7. The particle and wall dynamics are integrated 
using a Verlet algorithm. These discrete elements are 
coupled to a density matched fluid, described as a slowly 
varying continuum phase. The hydrodynamical fluctua- 



tions of the pores [T7] are neglected. The fluid velocity 
uf{y) and the fluid shear stress cr^ (y) profiles are deter- 
mined by averaging the equations governing the motion 
of the fluid over x and over time t, as proposed in [TS] 
(see Supplemental Material). 

The particles are submitted to four types of forces, (i) 
Upon contact, they interact with a viscoelastic force and 
with a Coulomb friction for relative tangential motion 
between particles at contact [T5H2T] . The model used 
for particle-particle interactions is identical to that pro- 
posed by Luding [21] . Quantities used in the model are 
expressed in terms of the grain density p, of the applied 
pressure P p , and of the mean grain diameter d. In this 
system of units, the normal spring constant is chosen suf- 
ficiently large (between 10 3 and 10 4 ) to reach the rigid 
regime in which the results do not depend on it. The 
Coulomb friction coefficient is chosen equal to fi p = 0.4, 
except for the inset of Fig. [2] which is obtained in the 
frictionless limit p p = 0. The other viscoelastic param- 
eters are chosen to lead to a restitution coefficient small 
enough (between 0.1 and 0.9) to get results that do not 
depend on it (see Supplementary Information) . (ii) They 
are submitted to a viscous drag force given by: 



(1) 



which involves the nonaffine particle velocity component, 
i.e. the fluid velocity minus the particle velocity 
u p , and where i is the particle label. This is based on 
the assumption that the particle based Reynolds num- 

p \u p — u? I c 



bers Re 



remains small, (iii) When the 



fluid presents a stress gradient, it exerts a resultant 
Archimedes force on the particle, which reads /" chl = 
4>(1 - 0) _1 /f' as (see Supplemental Material), (iv) Fi- 
nally, when particles are separated by a lubrication film, 
we include the extra stress as an interparticle force me- 
diated by the fluid [22]: 
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where hij is the gap between the particles labelled i and 
j, d^ = -j^j: is the effective grain diameter, mj and ijj 
are the normal and tangential unit vectors between the 
grains. S is a regularization length, chosen equal to 5% 
of particle diameter. In real suspensions, it can be either 
related to the slip length, to the grain roughness or to the 
scale over which grains are elastically deformed [55 . This 
lubrication interaction is truncated for hij > (di +dj)/A. 

As the fluid is described as a continuum phase in a 
steady state, inertial effects and nonaffine effects are en- 
tirely ascribed to the particle phase. This means that the 
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FIG. 3: Simulated data, a): <f> as a function of K (inset as a function of J), b): (i as a function of K (inset as a function of 
I 2 ), c): f T as a function of <f> (inset shows f p ). All three Figures are for various I / J with the same color coding as in Figure 
M Solid lines corresponds to fits according to Eqs. 6,7, f T = fJ./K and f p = l/K. Dashed lines show the critical values. 



density p only appears in the equation of motion for the 
grains and includes the added-mass effect. 

As obtained for dense granular flows [HI], the simu- 
lation is insensitive to microscopic parameters provided 
that the grains are hard enough. The state of the system 
is then characterised by the two dimensionless numbers 
I and J. In the following, we will rather use the Stokes 
number I 2 / J = pjd 2 /rjf and the rescaled confining pres- 
sure I/J= \/pP p d/rif. 

Transition from viscous to inertial regime. - Fig. [I] 
presents simulation results obtained at the same vol- 
ume fraction tf> ~ 0.78 by varying the rescaled confin- 
ing pressure I/J, where / is typically varied from 10 -3 ' 5 
to 10 -0 5 . It compares the contributions to the dissi- 
pated power of the different forces acting on the bulk of 
the suspension. This dissipation is balanced by the en- 
ergy brought through the boundaries of the element of 
suspension considered. While the dissipation due to the 
drag force is dominant at small I / J, the dissipation in 
the contacts becomes dominant at large I / J and the sys- 
tem resembles a dry granular flow (within the inclusion 
of the added mass effect inside the density p). The sys- 
tem therefore presents a transition from a viscous to an 
inertial regime, controlled by the rescaled pressure. It 
can be seen that the dissipation in the fluid, both in the 
pores and in the lubrication films, gives a subdominant 
contribution and varies like the contribution due to the 
drag force. In the following, we will therefore focus on 
results obtained without the lubrication forces. 

Looking at the inset of Fig. [l] one observes that the 
friction coefficient p defined as the ratio of the particle 
shear stress t p and confining pressure P p remains con- 
stant across the transition. This means that, at fixed <f>, 
the shear stress is controlled by pressure, with a multi- 
plicative factor insensitive to the nature of the dissipation 



mechanisms. Fig. [2] shows the friction coefficient p of the 
system as a function of the volume fraction 4> for different 
values of the number I /J. A good data collapse is ob- 
tained, when I/J is changed over five decades, showing 
that p is a sole function of <j>. Moreover, the data ob- 
tained for frictional (p p — 0.4) and frictionless (p p = 0) 
particles fall on the same master curve jT9| and differ 
only by the values of <f> c and p c . 

A single rheology across the transition. - It has been 
recently argued that, in the viscous quasistatic limit, tra- 
jectories are mostly controlled by geometric effects close 
to the jamming point, and do not depend much on the 
nature of the mechanisms dissipating energy |I0j . We hy- 
potheses here that the paths along which particles move 
do not vary much across the viscous/inertial crossover. In 
the viscous quasistatic limit, it was shown that nonaffine 
displacements control the enhanced dissipation close to 
jamming and take place over a time-scale vanishing as 
7 _1 (<p c — 0)' 3 ^ 2 , where j3 ~ 2 is the divergence exponent 
of the stresses [10l [11]. In the inertial regime, micro- 
scopic rearrangements take place over an inertial time 
scale d^p/PP HHHl], over which we assume energy is 
dissipated. Assuming further that these two time scales 
are proportional to each other, we find that the stresses 
must also diverge as (<f> c — in this regime. 

Under the assumption that, for a given volume fraction 
4>, the dissipation due to viscous effects and that due to 
grain binary interactions can simply be added, particle 
shear stress and confining pressure can be written as sums 
of linear (viscous) and quadratic (Bagnold) terms in 7 
[H 023 : 

r p = Ui^ivn + apd 2 ^ 2 ), (4) 
P p = f P W(Vfi + apd 2 ^ 2 ). (5) 

The best fit of p and 4>, functions of / and J, give a con- 
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stant value of a = 0.635 ± 0.009 [Figs.[3j(a),(b)]. Indeed, 
we obtain a collapse of all data when fi and <fr are plotted 
against K = J + ctl 2 ■ This supports our above hypoth- 
esises. Consistently, expressions Q and ^ give the two 
relations (j) = f^O-fK) and \i = f T { ( t > )/fp{4 1 )- Following 
empirical expressions proposed for </> and \i as functions 
of / or J in the cases of dry granular flows and dense sus- 
pension respectively p] HH |24] , we can generalize them 
using the number K as 



(6) 
(7) 



fi{K) = fj, c + 



<j>{K) =4> B - bVK, 



where <j> c = 0.8139 ± 0.0003 is the jamming volume frac- 
tion. The constants b, fi c , jip, and Kq are specific to the 
considered system. Here we find, b = 0.412 ± 0.006, /x c = 
0.277 ± 0.001, n F = 0.85 ± 0.01 and y/K^ = 0.29 ± 0.01. 
Combining the two constitutive laws we finally get 
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This expression is in good agreement with the data 
displayed in Fig. [2j Furthermore, as f T = fi/K and 
fp = 1/K, these two functions are predicted to diverge 
close to the jamming point as {4> c — 4>)~ 2 , as a consequence 
of Eq. This behavior is also very well supported by 
our data, as seen in Fig. [3^c). 

We have run simulations in which lubrication interac- 
tions between the grains are taken into account. They 
do not affect the qualitative results described above but 
slightly change the values of the constants. In particular, 
the exponent of the diverging behavior of both functions 
f T and f p is unchanged. This contradicts the claim of [9 
that the divergence would be in {<p c — (j))" 1 when lubri- 
cation forces are present. 

Discussion. - The above analysis shows a crossover 
from viscous to inertial flow at a Stokes number I 2 / J = 
jd 2 p s /r]f ~ 1/a. The suspension is therefore found to 
present shear thinning at controlled granular pressure 
or shear thickening at controlled volume fraction, when 
I 2 j J goes beyond this value. In the experiments of Boyer 
et al, the maximum value of the Stokes number can be 
estimated as 10~ 3 . This value is far below the inertial 
regime, and, consistently, all their rheological data col- 
lapse when using J as the single dimensionless parameter 
PQ. By contrast, Fall et al. report in their experiments 
a crossover between the two regimes, at a Stokes num- 
ber of 2 10 -3 [T2]. This value is three to four orders of 
magnitude lower than the predictions of our simulations. 
We hypothesize that this effect may result from nonlocal 
effects, as the base flow is heterogenous. The dominant 
influence of nonlocality has previously been observed in 
other heterogenous flows of dense suspensions , emul- 
sions [57] . and granular systems [2"5H3"U] . Our setup is 
insensitive to nonlocal effects as all studied quantities are 



homogenous in the central part of our shear cell. Nev- 
ertheless, many flows are heterogenous and it would be 
important to understand nonlocality in order to rational- 
ize even these. 

Newtonian fluids exhibit a transition from laminar to 
turbulent flow controlled by the Reynolds number based 
on the size of the flow and on the suspension viscosity. It 
is unlikely that dilute or even moderately concentrated 
suspensions would be an exception to this rule. As the 
jamming transition is approached {4> — > (f> c ), the suspen- 
sion viscosity diverges so that the Reynolds number van- 
ishes. The transition from the viscous to the inertial 
regime in dense suspension is thereby of a different nature 
than the transition from laminar to turbulent flow. In 
the former, both the Newtonian and Bagnoldian regimes 
are controlled by particle fluctuations with respect to 
the affine field. These fluctuations are controlled by the 
Stokes number, which is based on the grain diameter and 
the fluid viscosity rather than the suspension viscosity. 
Further studies are needed to investigate the transition 
from the inertial regime to the turbulent regime when 
the particle volume fraction is lowered. 

In this Letter, we have shown that the Newtonian rhe- 
ology of suspensions can be unified with the Bagnoldian 
shear-thickening regime for vanishing temperature. As 
pointed out recently by Ikeda et al. [31] . thermal and 
athermal suspensions seem physically distinct, making a 
unified description of glass and jamming transitions un- 
likely. Future studies will have to explain the difference 
in nature (if any) between mechanically induced fluctu- 
ations (i.e. nonaffine motion) at zero temperature and 
thermal fluctuations. 
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SUPPLEMENTAL MATERIAL) 
Solving the two-phase hydrodynamics model 

In the presence of particles occupying a volume fraction 
4>, the hydrodynamics is described by the two-phase flow 
Reynolds averaged Navier-Stokes equations [T]. For the 
fluid phase we have: 

P/(l-0^- = V-o-'-F + p/(l-0g. (9) 

Where pf denotes the density of the fluid phase, g the 
gravity field, cr^ the stress tensors for fluid phase (with 
crfj = —pfSij + rh ), F the coupling term between the 
fluid and particles, and where is the material derivate 
given by: 



Our model assumes a newtonian fluid (rl y = r/fdul/dy) 
in steady-state and without inertia. Assuming further 
that we are at zero gravity reduces Eq.(l) to: 

= V-er / -F. (11) 

The coupling term is given by the sum of the drag forces 
/^ rag and the Archimedes forces /^ rchl over all particles 
labelled k in a given volume dV: 

F = ^7 E (/f as +/r hi )- (12) 

kGdV 

The Archimedes force is given by /^ rchl = V P X/ -<r k , where 
V p is the particle volume and a[ is the fluid stress exerted 
on grain k. Applying Eqs. (??) and (??) over a single 
grain with dV ~V p /<fi one obtains: 

jjarchi J?drag /-, q\ 

tk "(W) /fc • { ' 

Which then leads to 

F = 1 L "V f dra s = 1 F dra s (141 

in agreement with Jackson [1] . Our formalism applies the 
two-phase formalism on a single grain scale and thereby 
accounts for some hydrodynamical fluctuations at this 
scale. 

The fluid velocity field is found by sampling the x- 
component of the coupling term and integrating twice: 

u£(y)=u£(0)+ f -( [ h {F x (y"))dy"+ri)dy' (15) 
Jo Vf v Jo ' 

Where (F x (y)) is space- and time-averaged over a thin 
horizontal region and 200 time-steps. A new u' x fluid 
profile is calculated every 200th time-step. To preserve 
no-slip boundary condition at the walls we added at each 
time a constant stress, Tq , such to ensure the no-slip 
boundary condition. The fluid profile was continuously 
iterated to convergence in the simulations using the two- 
phase coupling term as a feed-back mechanism for the 
fluid, with a damping mechanism which makes use of 
the average of the 10 4 last fluid profiles rather than the 
instantaneous fluid velocity profile. The numbers in our 
averaging protocol are chosen in such a way to have a 
good balance between velocity update/convergence and 
stability. We checked that other averaging numbers give 
the same results. 

Viscoelastic parameters 



(10) 



Units used in the model are expressed in terms of the 
grain density (p), granular pressure (P p ), and mean grain 



G 



diameter (d). With these scales k n = 10 3 , 10 4 (P p ) (gran- 
ular pressure), kt = 0.5 k n (tangential spring constant), 
P n « 1.33, 4.20 (^PPpd 2 ) (normal damping), and /3 t « 
0.94, 2.97 (y/PPpd 2 ) (tangential damping) for /i p = 0.4. 
This yields a coefficient of restitution e ~ 0.9. For fi p = 
runs with k n = 10 3 (P p ) and /3 n « 1.33, 23.43 (y/PPpd?) 
were also performed, corresponding to e ~ 0.9,0.1. Es- 
pecially at high shear rates (i.e. high I) for the friction- 
less grains one finds a dependence on fx vs J or <j) on the 
value of the coefficient of restitution e. We choose to only 
show data points which are independent of the choice of 
e (within the interval 0.1 to 0.9). This roughly corre- 



sponds to points with two or more contacts in average 
per grain in the frictionless case. We see this as a sig- 
nature of being in a liquid-like rather than in a gaseous 
regime. 
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